Alex Williams, 09/2014
This ipython notebook demonstrates how to use PyNeuron Toolbox, along with the NEURON python module, to visualize the propogation of an action potential down an axon and in a model with realistic morphology. NEURON allows you to make visualizations like this with "shape plots" (see this tutorial). This demo will show you how to use PyNeuron Toolbox to make similar plots directly in python.
First, we will run a short simulation of a cable with Hodgkin-Huxley currents. The details for running NEURON simulations in python can be found in this paper (Hines et al., 2009).
We will use two functions from PyNeuron Toolbox,ez_record and ez_convert, to quickly record the membrane potential in all compartments of our cell. Using these functions is completely optional to making shape plots. Briefly, ez_record records the membrane potential for all compartments, and ez_convert converts all the h.Vector() objects into numpy arrays. The result is a matrix v whose columns are the membrane potential traces for all compartments (the labels for each compartment are returned by record_v for reference). These functions can also be used to record other common state variables in Neuron, such as the intracellular calcium concentration.
# Unmyelinated cable with Hodgkin-Huxley channels
from __future__ import division
from neuron import h
import numpy as np
import pylab as plt
%matplotlib inline
# PyNeuronToolbox imports
from PyNeuronToolbox.record import ez_record, ez_convert
# create model
h.load_file('stdrun.hoc')
cell = h.Section()
cell.nseg = 11 # It is a good idea to have nseg be an odd number
cell.Ra = 35.4 # Ohm*cm
cell.insert('hh')
# create 3d structure
h.pt3dadd(0,0,0,1.0,sec=cell)
h.pt3dadd(1732,1732,1732,1.0,sec=cell)
# Specify current injection
stim = h.IClamp(1.0,sec=cell)
stim.delay = 1 # ms
stim.dur = 100 # ms
stim.amp = 0.2 # nA
# Segment positions, equall spaced from 0 to 1
seg_positions = np.linspace(0,1,cell.nseg)
# Use toolbox to record v
# ez_record records v in all compartments by default
(v,v_labels) = ez_record(h)
# Manually record time and current
t, I = h.Vector(), h.Vector()
t.record(h._ref_t)
I.record(stim._ref_i)
# Run the simulation
h.init()
h.tstop = 25
h.run()
# Use toolbox convert v into a numpy 2D array
v = ez_convert(v)
Next we import the function shapeplot from PyNeuronToolbox.morphology. Then we simply create a 3D axis in matplotlib and call the function. This creates a plot of our axon in three dimensions.
from PyNeuronToolbox.morphology import shapeplot
from mpl_toolkits.mplot3d import Axes3D
# Plot shapeplot
plt.figure(figsize=(6,6))
shapeax = plt.subplot(111, projection='3d')
shapeplot(h,shapeax,lw=4)
plt.show()
Now lets make it interactive to show the propagation of the action potential. PyNeuron Toolbox has a function, shapeplot_animate, which makes it very easy to do this. This function returns another function which can be used by matplotlib animation tools. For a overview of these tools, see this blog post. In this example, we will use JSAnimation which makes it easy to embed animations on a webpage, like this one! The resulting animation shows the action potential propogating downwards from the upper left.
Note: We are using membrane potential for this example, but the same functions could be used to visualize changes in intracellular calcium or other variables of interest.
# Animation imports
from JSAnimation import IPython_display
from matplotlib import animation
# PyNeuronToolbox function for easy animation
from PyNeuronToolbox.morphology import shapeplot_animate
# Plot shapeplot (as above)
fig = plt.figure(figsize=(6,6))
shapeax = plt.subplot(111, projection='3d')
shapelines = shapeplot(h,shapeax,lw=4)
# Make the animation
nframes = 100
shape_animate = shapeplot_animate(v,shapelines,nframes)
animation.FuncAnimation(fig, shape_animate, frames=nframes, interval=40, blit=True)
It is also very easy to extend this approach to make more sophisticated animations. The code below creates an animated shapeplot on the right, with membrane potential traces on the left.
fig = plt.figure(figsize=(14,10))
# Plot voltage trace for each segment
lines = []
for i in range(1,cell.nseg+1):
plt.subplot(cell.nseg+1,2,2*(cell.nseg+1-i)-1)
line, = plt.plot([], [], lw=2)
lines.append(line)
plt.ylim([-80,60])
plt.xlim([0,h.tstop])
plt.yticks(np.arange(-80,60,50))
plt.ylabel('v (mV)', fontweight='bold')
plt.xticks(np.arange(0,h.tstop,5))
plt.setp(plt.gca().get_xticklabels(), visible=False)
# Plot injected current in bottom trace
plt.subplot(cell.nseg+1,2,2*cell.nseg+1)
line, = plt.plot([], [], '-r', lw=2)
lines.append(line)
plt.xlim([0,h.tstop])
plt.ylim([-0.05,0.3])
plt.yticks([0, 0.2])
plt.xlabel('t (ms)', fontweight='bold')
plt.ylabel('i (nA)', fontweight='bold')
# Plot shapeplot
shapeax = plt.subplot(1,2,2, projection='3d')
shapelines = shapeplot(h,shapeax,lw=4)
# convert 't' and 'i' to lists
t = t.to_python()
I = I.to_python()
# Animate function -- reveal traces, and animate shape plot
def animate(i):
# Animate the shape plot
shape_animate(i)
# Add additional animations
i_t = int((i/nframes)*v.shape[0])
for s in range(v.shape[1]):
lines[s].set_data(t[0:i_t], v[0:i_t,s]) # voltage traces
lines[-1].set_data(t[0:i_t], I[0:i_t]) # current trace
# Make animation
nframes = 100
shape_animate = shapeplot_animate(v,shapelines,nframes)
animation.FuncAnimation(fig, animate, frames=nframes, interval=40, blit=True)